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Abstract 

A critique is presented of the frequently used Bruce- Wilding (BW) mixed-field scal- 
ing method for estimating the critical points of nonsymmetric model fluids from grand 
canonical simulation data. An explicit, systematic technique for implementing this 
method is set out thereby revealing clearly a fortunate, close cancelation of contribu- 
tions from the leading correction-to-scaling and thermal scaling functions that makes 
the method effective for Ising-type systems but which lacks a general theoretical basis. 
The BW approach does not allow for pressure mixing in the scaling fields which is es- 
sential for representing a Yang- Yang anomaly, namely, the divergence at criticality of 
the second temperature derivative, {d'^ fi^/dT'^), of the chemical potential fJ-a{T) on the 
phase boundary; but such behavior must be expected in realistic models. We show that 
allowance for pressure mixing does not alter the leading dependence of the critical tem- 
perature estimator, Tc(L), on the linear size, L, of the simulation box: this converges 
as L^^^^^y'^ when L ^ oo (where u ~ 0.6 and 9 ~ 0.5 are the correlation-length and 
leading correction critical exponents). On the other hand, the critical density estima- 
tor, Pc{L), gains a leading variation oc L"^^/*^ that dominates the previously claimed 
j^-(i-a)/u ^^gj.j^ (where a ~ 0.1 and f3 ~ 0.3 are the specific heat and coexistence 
curve exponents). Numerically, the BW method provides estimates of Tc consistent 
with those obtained from recently developed unbiased techniques that do not require 
prior knowledge of the universal order-parameter and energy distribution functions; 
however, BW estimates of the critical densities, pc, prove significantly less reliable. 
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I. Introduction and Overview 



The critical behavior of fluids, including the issues of exponents and universality class 
and the magnitudes of non-universal quantities, such as the location of critical points, has 
been much studied by Monte Carlo simulations of model systems with the aid of increasingly 
powerful computers. However, since simulations always deal with flnite systems, say of linear 
dimensions L, a precise calculation can do no more than reveal details of the rounding of the 
system's properties in the critical region: a sharp critical point is attained only in the limit 
L ^ oo. Thus the analysis of flnite-size effects by suitable scaling concepts and subsequent 
extrapolation is a key ingredient for the success of such investigations.^ 

A notable advance in this respect was made in 1992 when Bruce and Wilding^'^ (BW) 
developed a special finite-size mixeA-fi.eld scaling technique or "recipe'' that has been widely 
used in deriving critical parameter estimates for a variety of more-or-less realistic model 
fluids.^~^° In their approach the joint distribution function, li), of the (number) density 
p = N/V and the configurational energy density, m, for a finite-size system of volume V = L''' 
and dimensions L x L x ■ ■ ■ x L with periodic boundary conditions is calculated via grand 
canonical ensemble simulations of the model fluid at temperature T and chemical potential 
p. Within scaling theory BW then described this joint distribution of density and energy 
fluctuations in a flnite near-critical fluid by incorporating the mixing of two relevant scaling 
fields; in other words, the ordering fleld h (conjugate to the order parameter) and the thermal 
fleld t (conjugate to the critical energy density) were taken as linear combinations of the 
"bare" physical flelds, t ocT — Tc and h oc // — //c-^^~^^ This mixing of t and h into the scahng 
fields h and t is a crucial manifestation of the lack of particlc-hole symmetry. An analysis 
of a {d — 2)-dimensional Lennard- Jones fluid using Monte Carlo simulations conflrmed this 
mixing.^ The observed (or calculated) joint distribution was then related to the universal 
critical-point distribution for the Ising model, to which universality class "typical" fluids 
are generally believed to belong. It was found^ that the limiting critical behavior of the 
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density distribution assumed a scaling form, as expected; furthermore, the shape of the 
distribution could be quite well matched to the known universal order-parameter distribution 
function appropriate to the Ising class. This served to illuminate the underlying basis of 
the universality shared by simple fluids and Ising ferromagnets. By such a matching of 
the distribution functions of density and energy to the limiting universal functions — that 
were obtained a priori via computer simulations of simple, symmetric Ising models^'^^'^^ — 
Wilding and subsequent collaborators were able to estimate the critical parameters and the 
mixing coefficients for more general continuum model fluids. 

Although the BW method has proved convenient and straightforward and has been suc- 
cessfully applied to estimating the critical parameters in various models, it suffers from 
certain weaknesses. One of its weakest points is that various crucial features of the critical 
behavior of the model under the investigation must be known a priori. In particular, the 
method requires precise advance knowledge of the fixed-point distribution functions. In- 
deed, Camp and Patey^^ recently studied criticality in three-dimensional model fluids with 
algebraically decaying attractive pair interactions varying as —1/r^^'^ with exponent values 
cr = 3, 1, and 0.1, via grand canonical Monte Carlo simulations. When they applied the BW 
method, they found that the order-parameter distribution for cr = 3 could be well matched 
to the Ising flxed-point function; this is, in fact, consistent with the renormalization group 
theory conclusion that for cr = 3 the system belongs to the Ising universality class. ^^'^'^ 

On the other hand, they observed that for o" = 1 and 0.1 the order-parameter distributions 
could not be matched to the Ising fixed-point function so successfully. These observations 
might well have been anticipated since RG theory indicates that both these systems should 
behave classically, ^^'^^ i.e., display van der Waals critical exponents {a — 0, /3 = ^,7 = 
1, 1/ = |).^^ Camp and Patey then tried to match the computed distribution functions to the 
approximately known classical fixed-point function^^ but with no success. They attributed 
the failure of the matching to lack of a precisely known classical fixed-point function; but it 



3 



could well be due to the more singular corrections-to-scaling associated with slowly decaying 
power-law interactions or to the failure of a fortunate but seemingly accidental cancelation 
of such corrections found in Ising systems as explained below. Indeed, even if the fixed 
point function were accurately known, unresolved questions about the general validity and 
effectiveness of the Bruce- Wilding analysis would still remain, especially for fluids whose 
behavior deviates significantly from particle-hole symmetry such as polymer solutions and 
"Coulombic" vs. "solvophobic" electrolytes.^^ These questions will be addressed below. 

One significant and potentially serious aspect of the BW approach has come to the fore 
only recently. This pertains to the presence of a Yang- Yang anomaly^^ in the bulk critical 
behavior of the system of interest. To be explicit, suppose iia{T) is the chemical potential 
on the phase boundary at and below Tc : then a Yang- Yang anomaly is represented by the 
divergence (or, more generally, corresponding singular behavior) of the second derivative, 
{(P IJ,„ / dT'^) when T T^—. The strength and sign of a Yang- Yang anomaly is conveniently 
measured^°'^^ by TZ^, the limiting ratio of C,ji,{T) = —T{(fiJi^/dT'^) to the constant- volume 
specific heat, Cv{T; p — pc) (or, more generally, to its singular part). 

It must be emphasized that although Yang- Yang anomalies vanish identically by con- 
struction in simple fluid models with an underlying gas-liquid (or 'particle-hole') symmetry, 
such as the standard lattice-gases in which PaiT) is always an analytic function through 
Tc, they must appear in general, as Yang and Yang^^ argued nearly 40 years ago. Further- 
more, both experimental data^^ and unequivocal simulation evidence — for the hard-core 
squarc-wcll (HCSW) fluid"^'^'"^^ and the restricted primitive model (RPM) electrolyte^^ — 
point to significant nonvanishing magnitudes of the Yang- Yang ratios TZ^ (even though the 
experimental situation is not yet fully transparent^^). In addition, an exactly soluble class 
of "compressible cell gas" models explicitly exhibits Yang- Yang anomalies. 

But how do such anomalies impact the BW method? The answer is that BW specifically 
failed to allow for the possible mixing of the pressure p into the scaling fields h and i via a term 
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proportional to p — pc- It transpires, however, that this previously neglected feature^^'^^ of a 
"full" scaling theory is essentiaF*^'^^ to account for the presence of a Yang- Yang anomaly in 
the bulk critical behavior. Accordingly, one of the aims of the analysis presented below is to 
investigate how the introduction of pressure mixing, which entails two further, independent 
mixing coefficients, affects the precision and reliability of the BW method. 

To tackle these issues we will invoke a rather thorough analysis of the full bulk scal- 
ing theory that incorporates pressure mixing and is applicable to nonsymmetric systems. ^'^ 
Then we will appeal to a recent systematic extension of the theory to systems of finite-size.^^ 
However, even to critique the BW theory in the absence of pressure mixing, it has proved 
necessary to formulate the BW matching procedures more precisely. It seems, indeed, that 
neither concrete details of actual implementations of the BW approach nor an explicit anal- 
ysis of the relevant finite-size effects are available in the literature. Accordingly, in section 
III below, we revisit the BW approach (without pressure mixing) from first-principles and 
describe the essential details in a systematic and explicit manner. 

On this basis, but subject to at least one serious proviso, we confirm the principal 
claims. ^'^ In particular, the effective finite-size critical density, Pc{L), defined by BW ap- 
proaches Pc like 1/L^^~°''>/^ when L — > oo (where a and v are critical exponents for the 
specific heat and the correlation length^^). Likewise, the corresponding effective critical 
temperature, Tc{L), approaches Tc as l/L*^^"*"^^/^ (where 9 is the leading even correction-to- 
scaling exponent^^). Armed with this information, effective extrapolations to L = oo may 
be undertaken. In addition to verifying these convergence laws, we obtain the corresponding 
leading amplitudes of the deviations, ATc(L) and Apc{L), in terms of parameters that are 
known or given. 

We then ask: Does pressure mixing change these results and, if so, how? Can the 
approach be modified readily to allow for the changes? In answering these questions, we 
use the finite-size scaling form for the canonical free energy, f{p,T), derived in ref 37. 
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We conclude, in fact, that pressure mixing will generate a new term in Pc{L) varying as 
l/L'^f^/" (where /5 is the cocxistcncc-curve exponcnt^^) that dominates the previously derived 
1/L*^^~")/^ term for large L. For Ising type (d = 3)-dimensional fluids one has 2(5 /v ~ 1.03 
while [1 — a)/v 1.41 so the difference exponent is only about 0.38.^^ 

The amplitude of this L~'^^/'^ term is proportional to the coefficient j2 that specifies the 
degree to which the pressure mixes into the ordering field h: see eq 3 below. Thus, when 
pressure mixing is negligible, i.e., j2 is sufficiently small, the BW estimation procedure for 
(that assumes the asymptotic l/L*^^"")/*^ behavior) remains reasonable. However, one should 
expect the dominant 1 / L^^/'^ term to play a more significant role in highly asymmetric fluids 
like the RPM where pressure mixing seems likely to be larger. Furthermore, these two terms 
may compete which could well yield potentially misleading nonmonotonic behavior. In the 
case of T'c(L), our analysis indicates that pressure mixing does not alter the leading l/L^^+^^Z'^ 
term although its amplitude changes if j2 7^ 0; however, that does not invalidate the BW 
approach. 

At this point we should point out that to fulfill the need for bias-free methods of 
assessing fluid criticality and, in particular, to convincingly determine universality class 
without prior assumptions, novel approaches have recently been devised, analyzed, and 
implemented. ^^'^^'^^'^^ These techniques rest on the calculation of certain special finite-size 
loci in the {p,T) plane, notably the k-loci^'^ and, especially, the Q-loci,^^ pq{T;L), and the 
examination of various density fluctuation moments as functions of T and L on these loci. 
In this way the Ising character of criticality in the RPM electrolyte has been established, ^^'^^ 
and precise and, apparently, rather reliable estimates of Tc and pc have been found for the 
model. 

However, these new methods require more and more careful computation! Accordingly, 
for "every-day" or exploratory applications where Ising behavior may be reasonably pre- 
supposed, it is reasonable to ask what kind of error one should expect in applying the 
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convenient and relatively economical BW recipe to estimate the critical points of new model 
fluids. Indeed, the BW method does seem to give reliable estimates of when the univer- 
sality class of the model under investigation is of Ising character. Some concrete evidence 
for this conclusion is presented in Table 1, which gives estimates for the HCSW fluid and the 
RPM. One sees that for the rather symmetric HCSW fluid (with a very small Yang- Yang 
anomaly, TZ^ ~ —0.04) the central BW estimate for Tc is only same 1 part in 10^ higher 
than the unbiased value; indeed, the quoted uncertainties overlap. For the RPM (with a 
highly asymmetric coexistence curve and a signiflcant Yang- Yang anomaly, TZ^ ~ +0.26) 
the uncertainties no longer overlap; but the BW estimate is lower by less than 1 part in 10^. 
On the other hand, the critical densities provided by the BW method and their apparent 
precision are in much poorer agreement with the new unbiased estimates; for both models 
the BW values are higher, by 1% and 6%, respectively. One may suspect that a similar sit- 
uation will prevail more generally; but further precise simulations, for example, for polymer 
solution models, are needed to confirm such a verdict. 

The balance of this article is organized as follows: In section II we briefly recapitulate, 
for ease of reference and to deflne notation, the basic scaling theory with pressure mixing.^° 
The vital details of the BW method and concrete procedures for matching the distribution 
functions are formulated in section III. The limiting behavior of the basic estimators Pc{L) 
and Tc(L) is then discussed pointing out the fortunate accident (noted but not stressed 
by the original authors) that seems essential to its practical success. Section IV addresses 
the modifications of the BW analysis that allow for pressure mixing. In section V a brief 
summary and discussion is provided. 
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II. Full Scaling Theory for Fluids 



II. 1 Bulk thermodynamics 

Following refs 30 and 31 we consider a single-component fluid with a critical point at 
Pci i^c-i and Tc and with a critical density p^- Convenient dimensionless variables near the 
critical point are then 

p=ip-Pc)/PckBTc, p={p-po)/kBTc, t={T-T,)/T„ (1) 

where /cb is Boltzmann's constant. Neglecting terms beyond linear order^^ the nonlinear 
scaling fields, p, h, and i may be written 

P ^ p - kot - loll -\ , (2) 

h = p - ht - j2P -\ , (3) 

t = t-hfi- jip+---, (4) 

where, for future reference we mention that in the bulk or thermodynamic limit the choices 
in eq 1 leads to /q = 1 while is related to the critical point entropy. [See eq 3.22 of ref 30.] 
In these terms a general scaling hypothesis asserts that the thermodynamics near criticality 
can be described, at least asymptotically, by^^ 

p ^ QW-''w±{uh/\tf, uSW uS\''). (5) 

with ± corresponding to f ^ 0, where A = /3 + 7 is the gap exponent while a, j3 and 7 are 
the standard exponents. As mentioned above, 9 is the leading even correction-to-scaling 
exponent while 6^ is the corresponding leading odd exponent. The coefficients Q, U, U4 and 
U5 are nonuniversal amplitudes while, when suitably normalized, W± are two branches of 
a universal scaling function. The nonuniversal irrelevant scaling amplitudes, U4 and U5, will, 
in general, depend smoothly on the variables t, p, fi ox p, h and i?^'^^ Note that the BW 
assumption corresponds to setting ji = ^2 = (which does not, in and of itself, relate to 
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any special symmetry) ; thus their analysis requires the determination of only the two mixing 

coefficients ki and li. 

In line with the original BW analysis, it is appropriate to define two relevant critical 
densities or "operators," A4 and £, conjugate to the ordering field h and the thermal scahng 
field t in the underlying Hamiltonian. In the context of an Ising ferromagnet these may 
be identified as the fiuctuating magnetization and the fiuctuating energy density (or, more 
precisely, as the deviations from the respective mean values at criticality) . Thus we may 
obtain the mean values of these basic scaling densities via 

{SM) = {M - M.) = (^^^ ^, (5£) = (£-£e) = -(^)^, (6) 

where the subscripts denote mean values evaluated at criticality. (Note that in the previous 
analysis the notation adopted was p — (SM.) and s — —{SE). We may also define the 
thermodynamic reduced number density p and energy density u via^^ 

These are related to the number density, p = N/V, and the total entropy, S, by 

. p . S , . 

p = Pc + Ap = — , M = Mc + Am = - , (8) 

Pc VpckB 

where, as previously, is the number of particles in the system and V the volume. 

Now the relevant critical densities, {SA4) and (SS), can be expressed as nonlinear com- 
binations of the reduced energy and number densities: after some algebra one finds^^ 

/AA//\ ~ {I- jiko)Ap- {h+jilo)Au 

^ K-{n+nk,Wp+{n+nh)Au' 

_ (l-j2/o)A7i-(fci+j2fco)Ap . . 

^ ^ ~ ii'-(j2+Jifci)Ap+(ji+M)A^' ^ ' 

where the constant 

K = 1- kih - jiko - j2lo - j2koh - jikilo, (11) 
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reduces to 1 — kih in the absence of pressure mixing. Note that even though the relevant 
scahng fields are considered only up to linear order, the presence of the linear pressure- mixing 
makes the mean order-parameter {A4) and critical energy density (S) nonimear combinations 
of the number density p and the energy density it. This is in contrast to the linear forms 
assumed by BW for A4 and £ in terms of the fluctuating number and energy densities:^'^ 
however, in the absence of pressure mixing ( i.e., when ji — j'2 = 0) the two mean densities 
reduce to 

Ap-hAu Au-kiAp 

in accord with the BW formulation. 

Even in the presence of pressure-mixing, one can expand eqs 9 and 10 near the critical 
point up to linear order to obtain 

{5M) ^ K-^[il-jiko)Ap-ih+j,lo)Au], (13) 
{SS) ^ K-'[{l-j2lo)Au-{k,+j2ko)Ap]. (14) 

These forms are consistent with the expressions given by BW except for the modification of 
the coefficients by the pressure-mixing coefficients ji and j2. 

II. 2 Finite-size scaling 

To discuss the analysis of simulation data we need an appropriate, "complete" finite- 
size scaling formulation. To that end we follow ref 37 and first introduce the dimensionless 
size-scaled temperature and ordering field variables 

wl = OffLV-, zl = auhL^'''. (15) 

and the leading even and odd correction variables 

ZLA^a^L'^/'' and = OgL-^^/^ (16) 
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Here v is the correlation length exponent while ag, ax, 04, and 05 are nonuniversal metrical 
factors (depending only weakly on p, /x, T and L) of dimensions conjugate to the powers of 
L appearing in the definitions of wl, etc. The basic scaling hypothesis eq 5 is then extended 
to 

PcP ~ L~'^Y{wl, zl; zl4, zl5)- (17) 

See eq 2.2 of ref 37.^^ Apart from possible L-dependences of the scaling field^^ p, h and f, 
which, however, have no effects to the orders that will be relevant here,^'' the definitions 
in eqs 2, 3, and 4 still hold. The scaling function Y{w, z; Z4, Z5) should be universal (when 
suitably normalized) but must depend on the boundary conditions and, indeed, on the shape 
of the system domain. Accordingly, for concreteness we will presuppose cubic simulation 
boxes of edge length L with periodic boundary conditions. Since in a finite system the 
grand canonical pressure must be related analytically to the other fields, the scaling function 
Y{xl,- ■ ■) should be analytic for small values of all its arguments.^^ 

The bulk hmit may now be obtained formally by setting L — l/\a£t\'^ and letting L — > 
00. The form of eq 5 is then recaptured provided the hyperscaling relation du = 2 — a 
is accepted. ^'^ One further finds Q — \a£\'^~°' / Pc, which is dimensionless, U — ax/jagl^, 
U4 = a^la^l^, and f/5 = asjagl^^, while W±{z; z^, 2:5) = F(±l, z\ z^, z^). 

Now BW focus on the joint distribution function, Pl{A4,S), of the critical densities A4 
and £ and suppose this has an appropriate universal scaling form. Without doubt the mean 
values, {SA4)l and {S)l, in a finite system will still be given by the previous thermodynamic 
first derivatives — see eq 6 — provided p is replaced by the finite-size expression eq 17 
for p{p, 11, T; L).^^ However, the mean values clearly do not define the distribution function 
Pi[}A,8). On the other hand, from the successive derivatives {d^'^^^p/ dh^dV^) one can 
extract, in a standard way, the corresponding moments {{5Ai)"^{5EY) l for all {m,n) as 
functions of h and i (or, equivalently, as functions of p, and T). Then, from the full set of 
moments one may, at least in principle, reconstruct the corresponding distribution function 



11 



that is desired. Furthermore, it is clear that all the scaling properties of the set of moments 
will be inherited by the distribution function. 

If we introduce the scaled number and energy fluctuation variables via 

XL^5ML^''' /au and = ^^L^'-^^/Vag, (18) 

it follows that we may conclude 

Pl{M, £) ^ NlV{xl, Vl, wl, zl, zl4, ^ls), (19) 

where the normalization constant Af^ can play no special role. Apart from notation and 
the neglect of the leading odd correction variable zl5, this is precisely the finite-size scaling 
ansatz advanced by Wilding.^ (Of course we ourselves have, here and above, neglected all 
higher order correction variables ZLk for k > 5.23,31,37^ j-^. jg ^^^^ clear that the scaling function 

V{xl, ■ ■ ■) should be universal. At bulk criticality, where h = t = 0, so that wl = zl = 0, 
we can then write 

Pl{M, £) « NLV,{5ML^/''/aM, 5£L^^-"^'^ /a£- a^/L"\ a^jL'^'^) (20) 

where Vc{xl, ■ ■ •) embodies the universal, statistically scale- invariant critical density fluctu- 
ations characteristic of a specific universality class but also incorporates the leading inverse 
powers of L that must play a significant role in systems that are not very large. 



III. Analysis of the Bruce- Wilding Method 

The aim of the BW method is to match the observed fluctuation data for the particle 
number density, 5p, and energy density, 5u, to the expected near-critical ordering density 
distribution, that follows from 

PLMiM) = f PL{M,S)dS, (21) 
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on the basis of the scahng ansatz eq 20, and hkewise to the corresponding critical energy 
density distribution, pl^£{£). In essence the matching is to be accomphshed (i) by suitably 
choosing the two mixing parameters h and ki in the expressions 

5M oc 5p — li5u, 5£ (X 5u — ki5p, (22) 

(that underlie eq 12 for the mean values) and (ii) by determining L-dependent approxima- 
tions, Tc{L), Pc{L) and Uc{L), for the bulk values T^., Pc and Uc- Appropriate extrapolations 
on L should then yield optimal estimates for Tc and pc- 

The universal limiting (L — > oo) critical distributions, which we will callp^(x) andpg(y), 
where x and y are defined as in eq 18 (but with the subscripts L dropped for brevity) are 
assumed known a priori. Indeed, fluid-magnet universality implies that a critical fluid order 
parameter distribution, Pl^x(AI), should, when L oo, precisely match the universal func- 
tion p*j^{x) appropriate to the Ising universality class; but this can be found independently 
by careful Monte Carlo simulations of Ising models.^'^^'^^ Figure 1 shows this distribution as 
obtained by Wilding and Miiller^ via simulations of the {d — 3) -dimensional Ising model. The 
distribution function has been normalized to unit integrated weight while the nonuniversal 
amplitude qm is chosen so that the distribution has unit variance. Evidently Pm{x) has two 
symmetrical peaks at, say, x = ±x* of equal height where this value and x* should 

be universal: see Fig. 1. The same approach applies to the energy distribution Pl^£{S)- this 
should match p^y) which may also be found numerically:^ see Fig. 2. This distribution, 
which is unimodal but asymmetric, has been normalized similarly by choice of as. Again 
the location of the peak at, say, y — y*, and its height, p^^, must be universal: see Fig. 2. 
Other universahty classes, such as, e.g., the XY or Heisenberg classes, etc. may be expected 
to have broadly similar distributions but with distinct values of x*, y*, etc. 

Now let us first consider the BW theory, in the absence of pressure mixing, and formulate 
their matching procedures in an explicit manner. 
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III.l Matching conditions 

The finite-size ordering distribution, can be obtained computationally by 

transforming the joint distribution pl{p,u) as observed in a simulation (where here we re- 
gard p and u as the fluctuating finite-size values of density and energy), by using eqs 21 and 
22, while simultaneously transforming the chemical potential p and the temperature T by 
introducing the (as yet unknown) mixing coefficient h. As indicated, when L ^ oo, the 
resulting function, Pl,a4(A^), should match the universal distribution p%i{x) in terms of the 
variable 

X = Lf'I^M - M^)/aM- (23) 
Using the (known) critical exponents (5 and v for the universality class in question, matching 
Vlm^-^) to universal function pX^(x) at some fixed L is to be implemented by choosing 
five fitting parameters, which, following BW, we take as the estimators Tf.{L), PciL), h{L), 
aM{L) and M.c{L) since, when L — > oo, these should approach the corresponding bulk 
values. 

Wilding has given no details of how he actually performs the matching that he implements: 
this leaves some ambiguity when it comes to estimating finite-size corrections and studying 
the effects of pressure mixing. Accordingly, for precision we consider the following formal 
procedure: (a) At finite L and fixed T < Tc [and for values of p not too far from PaiT)] 
the measured distribution Pl,a^(A^) will normally be two-peaked. By adjusting p to, say, 
P^^\L) the peaks can be made of equal heights; But, in general they will not be symmet- 
rically disposed with respect to the intervening minimum. It should be possible to achieve 
this by (b) adjusting li and, generally, re-adjusting p to obtain values p^'^\L) and t{\L); 

{'2) 

(bl) the position of the minimum then identifies an estimate Aic {L) and (b2) one can 
then satisfy the universal peak-placement condition 

x^ = a-^{L) L^l^iMt^ - M,) = ±x\ (24) 

for the peaks at M.)^ {L) by adjusting aM to, say, aj^ (L): At this point, in general, the peak 

14 



heights will not satisfy the universal relation desired, namely, 

Pt ^ PLMiM^I^) = PMi±^*) ^ Pm^- (25) 
Nor will the height of the minimum satisfy the corresponding relation 

pI = PL,A.(-Mf ) = PMiO) = P*r^.- (26) 

See Fig. 1. (c) To satisfy these two conditions it should (ignoring finite-size and pressure 
mixing corrections) suffice to vary the temperature T. Changing T should yield a series 
of further estimates fj,^^\L), fj,^^\L), ■ ■ •, for fJ,c{L) and, likewise, for the other parameters. 
Unless precision is high, both eqs 25 and 26 should yield the same set of five fitting param- 
eters Tc(L), HciL), h{L), aM{L), and M.c{L): significant differences, if they arise, could be 
a consequence of either finite-size or pressure-mixing effects (or the 'known' universal dis- 
tribution could be in error). However, practical experience for Ising-type systems suggests 
that satisfactory fits within the numerical uncertainties can be obtained fairly readily; but, 
even this might change if data of significantly higher precision became available. 

To determine the remaining fitting parameters for the energy, namely, Sc{L), as{L), and 
ki{L), one then considers the variable 

y = L^'-»y^{S-S,)/a£, (27) 

in order to match the data for the finite-size energy distribution, Pl,£{^), to the known 
fixed-point function p£{y) of character as shown in Fig. 2. 

We also remark that if pressure mixing is absent, i.e., ji — j2 = 0, the fitting values Pc{L) 
and Uc{L) may also be determined by using eq 12. Notice then that the mixing parameters, 
Iq{L), and kQ{L) are simply equal to Pc{L) and —Uc{L), respectively, as follows from the 
scahng analysis.^-*^'^^ 



III. 2 Z-dependence of estimators for and fi^ 

To elucidate the L-dependence of the estimators Tc{L), Pc{L), etc., determined from the 
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fits at fixed L as set out above, let us integrate the joint distribution eq 19 over the energy 
fluctuations to obtain the order distribution, We may expand this distribution 

about the fixed point function (x) as 

Pl,m{M) « NM\p*Mi.^) + a£L'/''tp;{x) + aiL-"''pl{x) 

+ aML'^/'fiplix) + a,L-''/''pl{x) + ■■■], (28) 

where x is defined in eq 23 and the normalizating factor will play no role. The derivative 
scaling functions pl{x), pl{x), etc., should be universal. Note also that by the symmetry of 
the Ising (or similar) fixed point, the functions p%i{x), Pt{x) and pl{x) are symmetric under 
X —X, while pJX'x) and pl{x), which describe the contributions from the ordering field and 
the odd corrections-to-scaling, respectively, must be antisymmetric. 

The original BW approach assumed tacitly that throughout the matching procedure 
/i = is maintained; it also ignored the odd correction-to-scaling contribution. In this case, 
the order distribution Pl,x(A^) becomes symmetric about M. — M.c- This will, in fact, be 
true for systems displaying an Ising- type symmetry. However, for asymmetric (e.g., fiuid) 
systems, the ordering field h at bulk T^. and Pc will depend on the system size L. In fact, 
one can show^^ that in this case the ordering field h decays as 

h = aft/L(^-°+^)/'^ + • • • , (29) 

where, in the absence of pressure mixing, the amplitude dh is proportional to the mixing 
coefficient h. Note also that h generates a singular \t\^~°' term in the coexistence curve 
diameter. ^^'^^ Miiller and Wilding* also noticed this point in their study of an asymmetric 
binary polymer mixture and observed that the chemical potential deviation on the phase 
boundary decays with the same exponent {1 — a + ■y)/i'. The consequent presence of an 
antisymmetric contribution to Pl,m(A1), which varies as L^/^h oc L~'^^~"~^'>/'' then makes 
it more difficult to match the order distribution to the symmetric fixed point function for L 
not so large. This effect might cause some difficulties in applying the BW method to highly 
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asymmetric fluids even if pressure mixing may be neglected. 

Nevertheless, in order to elucidate the original BW method, let us initially ignore the 
contributions in eq 28 proportional to pl_{x) and pl{x) that arise from the nonzero value of 
h and the odd corrections-to-scahng. In that case, the order distribution Pl,m(-^) should 
be symmetric in A4 about the critical value /Ac- But then the rather precise collapse of 
PL,Mi-^) o^^o the flxed point function that has been achieved in practical simulations 

seems to require the effective cancelation of the contributions from the flnite non-zero value 
of the scaling field i by the leading even correction-to-scaling term: see the second and third 
terms in eq 28. Indeed, as noted by Wilding, Nicolaides and Bruce^'^^ the functional forms of 
the universal functions and pl{x) do seem to be similar in the Ising universality class, 

so that such a cancelation is feasible! However, this seems to be no more than a fortunate 
accident with only limited numerical support; no firm evidence as to why this should be true 
or should hold for other universality classes has been offered. Nevertheless, if one accepts this 
observation as a reasonable approximation, the cancelation of the two contributions yields 



as 



L^/"i + a4L-^/"7e* ~ 0, (30) 



where 71* is the (approximately) x-independent ratio of pl{x)/pl{x). The condition h = 0, 
which was accepted in the original BW method, leads via the scaling field definition eq 3, to 
(1 — kit where we are, here, neglecting j2. Substituting and using eq 4 (with ji — 0) leads to 

''^"-^ = (1 - hk)ae ' ^^'^ 

where is the true critical temperature. The exponent {1 + O)/^ here confirms the argu- 
ments of Wilding, Nicolaides and Bruce. ^'^^ For a nonzero value of ki, it also follows that 
the chemical potential estimator scales in the same way as Tc(L). We thus obtain explicitly 

- /.r - -^p^:^L-(^-^^)/^ (32) 

[1 — Kiii)a£ 

where //^ is the bulk critical chemical potential. Finally, recognizing eq 29 for h{L) , turns out 
not to change the leading behavior in the temperature and chemical potential estimators, 
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since, at least for (d = 3)-dimensional Ising-type systems, we have (1 — a + 7) ~ 2.13 > 
(1 + ^) ~ 1.5. But it must be emphasized again that this conclusion relies upon the surely 
approximate proportionality of pl{x) to pl{x). 

III. 3 Convergence of the energy and density estimators 

To study estimators for the true critical density, pc (= p^), and energy density, Uc 
(= M^), we consider the following equations derived by inverting eq 12, namely, 

Ap = (SM) + k {SS) , All = {d£) + h {SM). (33) 

In order to obtain Pc{L) and Uc{L) we must thus first determine the behavior of and 
£c{L)- Wilding and Miiller^'''' argued that the critical value of the ordering operator, namely, 
}Ac{L)., should not depend on the system size L, since the asymptotic order distribution 
Plm^-^) ^® symmetric in — {Ai) which fixes Aic- However, this cannot be strictly 
true if one recognizes the contributions from the nonzero value of h and the leading odd 
correction-to-scaling term. As mentioned, these contributions, proportional to p\{x) and 
p%{x), respectively, are antisymmetric in x. Therefore, for any finite L, a perfect collapse 
oi pl^m{M.) onto the symmetric, fixed point distribution p*j^{x) is not, in general, possible. 
Nevertheless, let us suppose that an optimal collapse of data has been achieved (say, by some 
least-squares procedure) yielding a best fit for the various estimators. How does M.c{L) then 
depend on LI 

The near-critical order distribution Pl,m{-^) have a local minimum &X, M. = M.c{L) 
in accord with the matching requirements. In the absence of p^(x) &ji<\p\{x) (and any higher 
order odd terms), this value should coincide with the limiting value /A^ (= M.'^) — which 
merely says that Pl,x(A4) has a local minimum at x = 0. However, the antisymmetric 
corrections, p^ix) and p%{x) must shift M.c{L) away from M.c so that Pl,m{^) '^i^^ ^^^^ ^ 
local minimum at some a; 7^ 0. To find this location, we expand the scaling functions in eq 



18 



28 about X — 0: according to their symmetries one may write 

PMi^) = + (34) 

and similarly for pX{x), while pX{x) and p%{x) generate only odd powers of x. At the mini- 
mum, the derivative of pl^x(AI), which takes the form 

+ 2aiL-^/''{pl^x + 2pl^x^ + • • •) 

+ aM~ahL-^'-''-^^'''{pl, + 3^3x2 + ■ ■ •) 

+ a,L-''l''{pl, + Sp^gx^ + •••) + •••, (35) 

must vanish. (Note that we have used eq 29.) On the other hand, the magnitude matching 
relation eq 26 yields the further condition 

= pI^)^^ + p(4)^4 ^ . . . ^ asL''%p;o + Pt2^^ + • • •) 
+ aiL-'^'''{plQ+pl^x^ + ■■■) 
+ aM~ahL-^'-''-^^'''{pl^x+pl^x^ + • • •) 

+ a^L-''^/''{pl^x + p^gx^ + ...) + .... (36) 

Solving these two conditions simultaneously for x and t yields 

i = _(a,pyp*„a,)/La+^)/'^ + ..., (37) 
X = -(a^a,p;;i/2pl2))/L(i-"-^)/'^ + --- 

-(a2Py2pf)/L^«/'^ + --.. (38) 

We may then note from the definition of TV in eq 30, that the p* ratio in the amplitude for 
i is simply pIq/pIq — TZ*. Prom eqs 23 and 38, we obtain 

5M{L) = MciL)-M^, 

« Ai/L^i-")/'^ + • • • + As/^^^^^'^^", (39) 
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where Ai — —a'j^dhPhi/'2pi' and A5 — —aM^bPti/'^P* - In contrast to the arguments of 
Wilding and Miiller,^'^ we thus find that Aic{L) does depend on L with a leading exponent 
— (1 — a)/v that, for Ising-type systems, is approximately —lAl?^ 

To obtain the L-dependence of the energy estimator £c{L)-, we first recall that the energy 
distribution p£;(y) is asymmetric with a maximum at y = y* 7^ 0: see Fig. 2. The matching 
of Pl,s{£) to pl{y) then yields 

5E{L) = £,{L) - ^ a£y7L(^-")/^ (40) 

Finally, from eq 33 we obtain 

p.{L)-p^ ^ {A, + hagy*)/L^'-^y'', (41) 
UoiL)-ur ^ {asy* + hA,)/L^'-"y\ (42) 

Notice that the leading exponent {1 — a)/^ agrees with the assertions of refs 6 and 7 so that 
the dependence of M.c{L) (contrary to the related claims^'^) does not actually disturb the 
anticipated asymptotic behavior. 

This completes our analysis of the BW approach when pressure mixing may be neglected. 

IV. Inclusion of Pressure Mixing 

As discussed in the Introduction, it is important in studying asymmetric fiuid criticality 
to understand and, if possible, to allow for the effects of pressure mixing in the BW ap- 
proach. In ref 37 a scaling expression for the bulk canonical free energy density, f{p,T), 
that incorporates pressure mixing is derived. On this basis a finite-size scaling form for the 
singular part of the canonical free energy density f{p, T; L) was advanced. A crucial feature 
is that pressure mixing introduces an extra correction term that vanishes as j2/L^^'^ when 
L — > 00. This contribution turns out to be antisymmetric in the ordering operator SA4. 
On noticing that Ib.\pl,m{-M)] becomes closely related to f{p,T) when L — > 00, we see that 
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the scaling ansatz for the order distribution Pl,m(-^) postulated by BW^'^ should better be 
modified to read, in expanded form [compare with eq 28] 

PlMM) ~ Mm [p*m{x) + aeL'/^ipK^) + j2a,L-f"''p]{x) + a^L-'/^'plix) 

+ aML^'^'hplix) + a,L-'''^pl{x) + •••], (43) 

where the new, i.e., the third term entails the scaling function p*{x) which should be uni- 
versal and antisymmetric in x while the amplitude, a^, is nonuniversal. In fact, this new 
contribution dominates all subsequent corrections when the exponents take Ising or similar 
values (where we appeal to eq 29 to see that L^/^h oc L~^^~'^~^^/^). 

The presence of this pressure mixing term evidently raises a further question about the 
validity of the BW method. When j2 is small, one may well still obtain good matching of 
the observed distribution PL,Mi-^) symmetric fixed point function p^(x) even for 

relatively small system sizes. The 2D Lennard-Jones fiuid may represent such a case, since 
Wilding^ observed that Pl,ai(A4) could be well symmetrized and readily matched to p*j^{x). 
However, if j'2 is sufficiently large, one should not ignore its contribution: then symmetriza- 
tion of Pl.m{-M.) should be feasible only in an approximate way even for relatively large 
L. Indeed, Caillol, Levesque and Weis^^ performed Monte Carlo simulations on equicharged 
hard-spheres (i.e., the RPM electrolyte) and observed that their data for pi^_v((A^) could be 
matched to the Ising distribution p^(a;) only for large L; for small L they were unable to 
symmetrize via the BW procedure. They attributed this failure to poor data sampling in 
the low density region of their smaller systems; but it would seem that significant pressure 
mixing in the modeP^ could well be the primary cause of the observed asymmetry although 
the antisymmetric contribution due to h, which varies hke L~'^^~°'~^^/^ [see eq 29], may also 
be a factor. Further simulations to resolve these possibilities would, thus, be interesting. 
Nevertheless, P%^m{-^) should always asymptotically approach p%[{x) when L — > 00. In 
practice therefore, one may still be able to match Pl,>i(A1) to P*m{x) within tolerable pre- 
cision for large enough L and thence derive best-fit estimators via the BW recipe. How will 
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these then depend on the system size LI 

Before addressing this question, we will revisit the matching conditions described in 
section III.l. As demonstrated by eqs 9 and 10, we must expect the critical ordering and 
energy densities, 6M. and 5£, to actually be non/mear combinations of the density and energy 
fluctuations, 5p and 5u. Near enough to the critical point, when the typical deviations are 
small, however, linear BW relations, following eqs 13 and 14, should still become valid. But 
we expect such linear relations to be inadequate further from criticality. 

Even if one can reach regimes where the linear relations are valid, however, the matching 
procedure should be more complicated, in order to accommodate the two extra unknown 
parameters, ji and j2- Here we propose an approach which, as far as possible, adopts the 
steps presented previously in section III.l: (i) First suppose ji — j2 — ^ and proceed through 
steps (a)-(c) in order to match the data for Pl,jw(A1) as well as feasible to P*m{x) and so 
obtain the first round of estimators Tc^\L), jj,i^\L), a^^{L), and Aii^\L). The 

remaining parameters, Sc^\l), k^\L) and ot£^(-^) can be obtained similarly by matching 
the energy operator distribution pl^£{E) to the fixed point distribution p^d/) as explained 
before. We should also recall the relations lo^^{L) = pi\L) and k^\L) = —u^\l). At 
this stage, however, one may still observe differences between the mixed data distributions 
and the fixed point distributions, especially further from criticality. (ii) Knowing via eq 43 
that j2 relates to asymmetry in Pl,.m(A^), we now suppose, as a tentative approximation, 
that the fluctuating critical densities, 5A4 and 5£, are related to the observable fluctuations, 
6p and 6u, via nonlinear relations that parallel eqs 9 and 10 (i.e., obtained by removing 
the expectation brackets and replacing the A's by S^s). Then, on first retaining the setting 
ji — 0, one can attempt to adjust j2 to obtain a value, say j^\L), that provides a better 
match oi Pl,m{-M.) to Pj^{x). Next (iii) one may adjust ji in the nonlinear relations to 
improve the matching of Pl,£{^) to p^iy): this should yield a value j^\L). (iv) Now one 
can set ji — j^\L) and j2 — j^\L) and recalculate the order distribution pL,Mi-^)' 
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is likely to observe some new discrepancies near the local minimum and the two maxima. 
Accordingly, one can return to step (i) but now with fixed ji = j[^\L) and j2 = j2^\L), and 
iterate the procedure to find a new set of parameters, say Tc'^\l), etc. On repeating 

these steps, one should be led to stable values for all the parameters. Nevertheless, as a 
practical matter one may reasonably question the robustness of this procedure (which we 
have not ourselves attempted to implement) . 

To obtain an assessment of the effect of pressure mixing on the convergence of the BW 
procedure, however, it suffices to suppose that we have achieved a good matching for the 
distribution functions pz,,x(Al) and pL,f (^) and have in hand a satisfactory fit of the critical 
parameters. To understand the L-dependence, we first expand the new universal function 
p*{x) in eq 43 as 

p*{x) ^ p*^x + p%x^ + ■ ■ ■ , (44) 

while the other functions appearing in eq 43 may be expanded just as in Sec. III. 3: see eq 
34. The local minimum of Pl,m(-^) "lust then satisfy 

= 2pf^x + 4pl^)x' + • • • + ^asL^'Hipl^x + 2p*^x^ + • • •) 

+ j2%^"^/"(p*i + 3p*3x2 + •••) + 2aiL-'/''{pl^x + 2pl,x'' + • • •) 

+ a^a,L-(^-"-^)/'^(p*, + 3p*3X^ + •••) + a,L-'^/^(pl, + 3pl,x' + ■■■) + ■■■. (45) 

while the matching condition, PL,Mi-^c) — P^^, yields a further condition corresponding to 

eq 36, namely, 

= +pWx^ + ... + asL'inipl^ + p^^x^ + ■ ■ ■) 

+ i2%i.-^/"(p*iX + p*3X=' + •••) + a^L-^'^'iplQ + pI^x'' + • • •) 

+ aMahL-^^-''~^^'''{pl^x + pl^x"" + •••) + a^L-'^/^pl.x + pl,x^ + ...) + .... (46) 

Solving these two equations simultaneously for x and t yields 

2p(2)a; = - j2a,p*jL^^^ - aM~aHPlJ L^'^''-^^''' - 2j2a,a,p*,plJ L^^+'^l^ + • ■ ■ , (47) 
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while, as regards the leading term, t is still given by eq 37. Note that the L"^^/^ term, 
that originally appeared in eq 38 (when ji = j2 = 0) still arises but only as a higher order 
correction not displayed here. [As previously, we assume tacitly that the exponent values lie 
within the normal range of 0{n) fixed points.] 

To derive the modified form for Tc[L), we first rewrite eq 2, the definition for p, in the 
form 

P^kot + lofi + p^ . (48) 

On using eq 3, the definition of similarly, the result in eq 29 then yields 

H = h±^lt + - + ,^^L-(^-"+-)/'^ + . . . , (49) 

1 - j2io 1 - J2h 1 - J2k 

which, on substitution, gives the leading order result 

P=f^o + ^o^^^lt + ---, (50) 



\ 1 - i2^o y 

from which we have dropped the terms varying as p oc L~^'^~°'^/^ (see ref 37) and L~^^~°''^'^'>^^ 
which enter only as higher order corrections. We may now substitute these results into eq 
37 and use eq 4 for t to find 

^c(^) = ^^^^^ - + . . . , (51) 

where the mixing coefficient combination is 

T = 1 - ji/co - {h + jilo)iki+j2ko)/il-j2lo), (52) 

while TZ* = p\q/pIq as in eq 38. The chemical potential may be obtained by substitution in 
eq 49 which yields 

(1 - j2lo)Ta£ 

Note that pressure mixing does not alter the leading exponent but does change the amplitude. 

Finally, the critical order estimator M.c{L) can be obtained from eq 47 which leads to 
the replacement of eq 39 by 

SM ^ -j2^2/^'^/" + Ai/L(i-")/'^ - j2^4/^^'^+^^/", (54) 
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where A2 — aMO-jP*jJ'^P* and — 204^43^2, while Ai is given after eq 39. On the other 
hand, to leading order, eq 40 remains valid for the energy estimator, 8c{L). To complete the 
calculation we now invert eqs 9 and 10 to obtain, up to linear order, 

Ap = {l-j2lo)5M + {h + 3^lo)5£, (55) 
= {l-hko)5E + {k^+j2kQ)5M. (56) 

Appealing to eqs 40 and 54 then yields our main conclusions, namely, 

Pc(L)-Pc(oo) = 5,L-2'3/'^ + ^pL-(i-")/- + B4pL-(2/5+')/'^ + ..., (57) 
u,{L) - tie(oo) = BuL-^P'^ + AuL-^^-'^^l^ + S4„L-(W)/- + . . . ^ (58) 

where the coefficients are 

Bp = -j2(l-j2/o)^2, Ap = {l-j2h)A^ + {li+jih)a£y\ (59) 
Bn = -i2(A:i+j2A;o)^2, A^^{1- jika)asy* + {h+j2h)Ai, (60) 
Bip 2aiP22-Bp and B^y, = 2aiP22-B«. (61) 

When the pressure-mixing coefficient j2 vanishes, the leading L~'^^/'^ terms drop out but 
the L~(^~")/'' terms remain; in that case one regains the original BW exponents for pc{L) 
and Uc{L) although the amplitudes now depend on ji. Evidently pressure mixing may 
significantly slow the rate of convergence in estimating p^. and u^- In practice it may be more 
significant that the exponents 2/3/i/ ~ 1.03, (1 — q;)/i/ ~ 1.41, and (2/3 + ^)/i/ ~ 1.86 (using 
Ising values) are fairly close in magnitude so that if the successive terms are of different 
sign and thus compete, reliable extrapolation may be seriously hampered. This could be the 
cause of the misleading BW error estimate for seen in Table 1 for the RPM. 

V. Conclusions 

In summary we have examined critically the Bruce- Wilding extrapolation method that, 
in the past, has been widely applied to various model systems since it usually provides a 
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straightforward and apparently reliable way of estimating critical parameters from finite- 
size data. We first analyzed in some detail the original BW method, that neglects pressure 
mixing in the scaling fields. Effective critical parameter estimators, Tc{L), Pc{L), etc., can be 
obtained by matching the numerically measured distribution functions Pl,jw(A4) and pl^s{^) 
to the fixed point functions {x) and {y) , respectively. We provided a precise specification 
for implementing the matching procedure that generates satisfactory fits and yields explicit 
values for the critical parameter estimators that can be investigated analytically. The finite- 
size behavior of the estimators, Tc{L), Pc{L), etc., was then derived. The principal BW 
claims, namely, that AT'c(-^) and AfidL) decay like l/L^^+^^Z'^ when L ^ oo while Apc{L) 
and Auc{L) vanish as Ij L^^~'^)/^ were confirmed with explicit expressions for the amplitudes: 
see eqs 31, 32, 41, and 42. 

When pressure mixing is allowed for, however, the ordering operator distribution, 
contains extra antisymmetric corrections that vanish only as l/L^^^; as a result these dom- 
inate all other corrections to scaling. Consequently, the matching of PL,Mi-^) sym- 
metric fixed point function (x) can be achieved only in an approximate way for finite L if 
one follows the BW recipe. An extension of the BW procedure that makes some allowance 
for pressure mixing was proposed but has not been tested. Nevertheless, by assuming that 
an acceptable matching can be realized, we demonstrated that pressure mixing does not alter 
the leading l/L^^'^^^^'^ term in the effective critical temperature, Tc(L): see eq 51 (although 
the amplitude is changed). On the other hand, the effective critical (number) density, pdL), 
and energy density, Uc{L), contain new, terms with amplitudes proportional to 

the pressure- mixing coefficient j2- Furthermore, these terms asymptotically dominate those 
previously identified: see eqs 57 and 58. 

In conclusion, the Bruce- Wilding method may still be regarded as a useful practical 
guide to the extrapolation of finite-size simulation data for systems that do not deviate far 
from symmetry and that may with confidence be expected to fall within the Ising critical 
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universality class: see Table 1 for some indications of its precision and numerical short com- 
ings. In other cases, however, the first issue of concern is that one lacks numerically reliable 
universal fixed-point distributions for the order parameter and energy that are essential for 
the method. Indeed, as outlined in the Introduction, Camp and Patey^^ implemented the 
BW method in a study of liquid-gas criticality in model fiuids with algebraically decaying 
attractive pair interactions, and encountered just this problem! Nevertheless, even if the 
required universal distributions were accurately known, one could not reasonably expect to 
benefit from the fortunate but apparently quite accidental cancelation of the thermal and 
leading even-correction scaling functions which have greatly assisted practical BW calcula- 
tions for Ising systems: see the discussion in association with eq 30 (where the ratio TZ* was 
introduced) . 

Beyond that, even when Ising criticality may be confidently assumed, the presence of 
both odd-order correction terms and significant pressure mixing, must be expected if, as in 
the case of the hard-sphere ionic fiuid models, the observed asymmetries are not small. 
In such cases, even when the BW matching recipe can be implemented satisfactorily, the 
results will very likely be distorted (relative to the more symmetric cases); consequently, 
the subsequent extrapolations must be regarded with increased caution and assessed as less 
reliable: see Table 1. 

While we have sketched one iterative BW-type method that could allow for pressure 
mixing, one might consider further elaborations of the BW approach — for example, by 
directly examining and "tuning out" the scaled cross correlations of the number density, p, 
and the configurational energy density, u. However, in the light of the recently developed 
bias-free procedures, involving the Q-loci and related estimators^^'^^'^^ (which seem to work 
rather reliably for at least some highly asymmetric systems and do not require detailed prior 
knowledge), attempts to further extend the BW approach do not seem warranted at present. 
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Table 1: Estimates for the reduced critical parameters of the hard-core square-well (HCSW) 
fluid and of the restricted primitive model (RPM) electrolyte via the Bruce- Wilding (BW) 
method^"^'^* and via bias-free methods. ^"^'^^ Parentheses denote stated uncertainties in the 
last decimal place quoted. 





HCSW 


RPM 


C 


Pi 


rp* 
C 


Pi 


BW 


1.2180(2) 


0.310(1) 


0.05065(2) 


0.084(1) 


Bias-free 


1.2179(3) 


0.3067(4) 


0.05069(2) 


0.0790(25) 



32 



Figure Captions 



Figure 1: The universal critical-point order-parameter distribution, as a function of 

the scaled order parameter x = aJ^L^^^'^dAd, for the Ising universality class, as calculated 
by Wilding and Miiller^ via Monte Carlo simulations of the {d — 3)-dimensional Ising model. 
The nonuniversal amplitude qm has been chosen so that the distribution has unit variance; 
the peaks of height — 0.4267 are then located at x = ±x* with x* ~ 1.1801, while the 
height of the minimum at x = is ~ 0.1904. 

Figure 2: The universal critical-point energy distribution, P£{y), for the Ising universality 
class as a function of the scaled energy y = ag^L^^^^'^^^SS, as calculated by Wilding and 
Miiller,® selecting the nonuniversal amplitude as so that the distribution has unit variance. 
The single peak of height p"^^ ~ 0.3981 occurs at y* ~ —0.3966. 
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